/*
 heat-spectral lib. solves time-dependent PDEs using SEM-NI method.
 Copyright (C) 2009, Altieri Michele and Stefanoni Andrea
 
 This program is free software: you can redistribute it and/or modify
 it under the terms of the GNU General Public License as published by
 the Free Software Foundation, either version 3 of the License, or
 (at your option) any later version.
 
 This program is distributed in the hope that it will be useful,
 but WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 GNU General Public License for more details.
 
 You should have received a copy of the GNU General Public License
 along with this program. If not, see <http://www.gnu.org/licenses/>.
*/

#ifndef ADV_HPP
#define ADV_HPP

#include<cmath>
#include<spectral>

// ******************************************************************
// Diffusion-Advection Problem
// ******************************************************************
// u_t - \mu \Delta u + b\cdot\Nabla u = 0	in \Omega \times [0,T]
//
// u(x,0) = u0					in \Omega
// \Nabla u(x,t) \cdot n = 0			on \Gamma_N
// ******************************************************************

#define ADV_PROBLEM AdvBilinear<spectral::Mesh>, Functional<spectral::Mesh,AdvFunc>


struct AdvFunc {
    double operator()(double x,double y) const {
	return 0;
    }
};

struct Advb1 {
    double operator()(double x, double y) const {
	return -y*std::sin(M_PI*std::sqrt(x*x+y*y));
    }
};

struct Advb2 {
    double operator()(double x, double y) const {
	return x*std::sin(M_PI*std::sqrt(x*x+y*y));
    }
};

template<class MESH>
class AdvBilinear {
    public:
	AdvBilinear(MESH const & me): me(me) {
	    for(int i=0;i<4;++i) {
		bcond[i]=spectral::neumann;
	    }
	}

	template<class U,class V> inline
	    double operator()(U const & u,V const & v,int n) const {
		return
		    + 0.001*me.prodL2(me.Dx(u,n),me.Dx(v,n),n)
		    + 0.001*me.prodL2(me.Dy(u,n),me.Dy(v,n),n)
		    + 0.1*me.prodL2(me.restriction(b1,n),me.Dx(u,n),v,n)
		    + 0.1*me.prodL2(me.restriction(b2,n),me.Dy(u,n),v,n);
	    }

	// Initial data
	struct initial{
	    double operator()(double const & x, double const & y) const {
		return std::pow(std::sin(2*M_PI*x)*std::sin(2*M_PI*y),6);
	    }
	};

	// Boundary conditions
	spectral::bctype bcond[4];

	// Function for the Dirichlet boundary condition
	struct dirichlet{
	    double operator()(double const & x, double const & y) const {
		return 0;
	    }
	};

	// Function for the Neumann boundary condition
	struct neumann{
	    double operator()(double const & x, double const & y) const {
		return 0;
	    }
	};

	// Methods on
	template<unsigned int N, typename T=void>
	struct on : public neumann {};

    private:
	MESH const & me;
	Advb1 b1;
	Advb2 b2;
};

const spectral::Point<double> LA(-2.0,-2.0);
const spectral::Point<double> LB(2.0,2.0);

#endif
